clear;clc;close all
addpath('./util');

%% 读取光片图像
LS = loadtiff('data/LS.tif');
%LS =loadtiff('G:\20251103\SP\LS\1\Image_1_51.tif');
%LS = apply_window(LS, 224, 803);
%figure;sliceViewer(LS,[]);
size(LS)

%下采样z，旋转，padding
LS = rot90(LS, 1);
LS = flip(LS, 2);

%% 读取光场图像
LF = loadtiff('data/LF.tif');
size(LF)

%% 对准xy

im_moving = squeeze(max(LS(:,:,20-10:20+10), [], 3));
im_fixed  = squeeze(max(LF(:,:,20), [], 3));

im_moving_adjust = imadjust(im_moving);
im_fixed_adjust = imadjust(im_fixed);

%手工粗配准
[moving_points,fixed_points] = cpselect(im_moving_adjust,im_fixed_adjust,'Wait',true);
t_concord = fitgeotrans(moving_points,fixed_points,'similarity');
fixed_ref = imref2d(size(im_fixed)); %relate intrinsic and world coordinates
moving_registered = imwarp(im_moving,t_concord,'OutputView',fixed_ref);

figure, imshowpair(imadjust(moving_registered),imadjust(im_fixed),'falsecolor'); title('人工配准结果');

%提高配准精度
[optimizer, metric] = imregconfig('multimodal');
optimizer.InitialRadius = optimizer.InitialRadius*2;
optimizer.MaximumIterations = 300;

movingRegisteredAdjustedInitialRadius = imregister(im_moving, im_fixed, 'affine', optimizer, metric, ...
    'InitialTransformation',t_concord);

% 找到xy的变换
tform_xy = imregtform(im_moving, im_fixed, 'affine', optimizer, metric, ...
    'InitialTransformation',t_concord);

figure, imshowpair(imadjust(movingRegisteredAdjustedInitialRadius), imadjust(im_fixed),'falsecolor'); title('精细配准结果');

% 保存XY配准结果
save('./run/M20_3Dtransfrom.mat', 'fixed_points','tform_xy');


%% 将LS，xy变换到LF

fixed_ref = imref2d([size(LF, 1), size(LF,2)]);
tLStiff_M20_beads = transform_LS2LF(LS, tform_xy, fixed_ref);


%% 方法一

%% 展示XY配准效果

im_moving = squeeze(max(tLStiff_M20_beads(:,:,:), [], 3));
im_fixed  = squeeze(max(LF(:,:,[1:30 70:101]), [], 3));
im_moving_adjust = imadjust(im_moving);
im_fixed_adjust = imadjust(im_fixed);

figure, imshowpair(im_moving_adjust, im_fixed_adjust,'falsecolor');  axis image;

pts = roi_pick_points(50);  %手动选取点进行Z轴配准
disp(pts);
%%
Ixz=imadjust(squeeze(max(tLStiff_M20_beads, [], 1)));
Iyz=imadjust(squeeze(max(tLStiff_M20_beads, [], 2)));
Ixz=Ixz';
Iyz=Iyz';

X_list=pts(:,1);
Y_list=pts(:,2);

z_list = pick_z_from_XZ_YZ_concat_rclick_batch(Ixz, Iyz, X_list, Y_list, 10, 'AutoInit', true);

disp(z_list);   % N×1，每个点对对应一个 Z

%%
LFxz=imadjust(squeeze(max(LF(:,:,[1:30 70:101]), [], 1)));
LFyz=imadjust(squeeze(max(LF(:,:,[1:30 70:101]), [], 2)));
LFxz=LFxz';
LFyz=LFyz';

LF_X_list=pts(:,1);
LF_Y_list=pts(:,2);

LF_z_list = pick_z_from_XZ_YZ_concat_rclick_batch(LFxz, LFyz, LF_X_list, LF_Y_list, 30, 'AutoInit', true);

disp(LF_z_list);   % N×1，每个点对对应一个 Z

%%

%对找到的beads的z轴坐标进行线性拟合
fit_xx=LF_z_list;
fit_yy=z_list;
p = polyfit(fit_xx,fit_yy,1); 
z_fit = polyval(p, fit_xx);

figure;
plot(fit_xx, fit_yy,'o' ,fit_xx ,z_fit ,'-');
legend('data','linear fit');



%%  方法二
% 对准xz

im_moving = squeeze(max(LS(:,:,:), [], 1));
im_fixed  = squeeze(max(LF(:,:,:), [], 1));

writetiff(im_moving,'./data/LSxz.tif');
writetiff(im_fixed,'./data/LFxz.tif');

im_moving_adjust = imadjust(im_moving);
im_fixed_adjust = imadjust(im_fixed);

%手工粗配准
[moving_points,fixed_points] = cpselect(im_moving_adjust,im_fixed_adjust,'Wait',true);
%t_concord = fitgeotrans(moving_points,fixed_points,'similarity');
%fixed_ref = imref2d(size(im_fixed)); %relate intrinsic and world coordinates
%moving_registered = imwarp(im_moving,t_concord,'OutputView',fixed_ref);

%figure, imshowpair(imadjust(moving_registered),imadjust(im_fixed),'falsecolor'); title('人工配准结果');